!(dirname(@__FILE__) in LOAD_PATH) && push!(LOAD_PATH, dirname(@__FILE__))
using DSP, Compat, Compat.Test, FilterTestHelpers

#
# Butterworth filter prototype
#

# Poles of 20 pole Butterworth filter prototype from MATLAB 2013b (buttap(20))
matlab_p = [-0.07845909572784487+0.996917333733128im,-0.07845909572784487-0.996917333733128im,-0.2334453638559053+0.9723699203976767im,-0.2334453638559053-0.9723699203976767im,-0.3826834323650897+0.9238795325112867im,-0.3826834323650897-0.9238795325112867im,-0.5224985647159488+0.8526401643540923im,-0.5224985647159488-0.8526401643540923im,-0.6494480483301835+0.760405965600031im,-0.6494480483301835-0.760405965600031im,-0.7604059656000308+0.6494480483301838im,-0.7604059656000308-0.6494480483301838im,-0.8526401643540922+0.5224985647159489im,-0.8526401643540922-0.5224985647159489im,-0.9238795325112867+0.3826834323650899im,-0.9238795325112867-0.3826834323650899im,-0.9723699203976766+0.2334453638559055im,-0.9723699203976766-0.2334453638559055im,-0.996917333733128+0.07845909572784507im,-0.996917333733128-0.07845909572784507im]
matlab_butter = ZeroPoleGain(Union{}[], matlab_p, 1)

# Test that our answers are close to MATLAB's and at least as accurate
butter = Butterworth(20)
zpkfilter_eq(butter, matlab_butter)
zpkfilter_accuracy(butter, matlab_butter, Butterworth(BigFloat, 20))

# Poles of 19 pole Butterworth filter prototype from MATLAB 2013b (buttap(19))
matlab_p = [-0.08257934547233227+0.9965844930066698im,-0.08257934547233227-0.9965844930066698im,-0.2454854871407991+0.9694002659393304im,-0.2454854871407991-0.9694002659393304im,-0.4016954246529694+0.9157733266550574im,-0.4016954246529694-0.9157733266550574im,-0.5469481581224267+0.8371664782625287im,-0.5469481581224267-0.8371664782625287im,-0.6772815716257409+0.7357239106731317im,-0.6772815716257409-0.7357239106731317im,-0.7891405093963935+0.6142127126896679im,-0.7891405093963935-0.6142127126896679im,-0.8794737512064892+0.4759473930370733im,-0.8794737512064892-0.4759473930370733im,-0.9458172417006346+0.3246994692046836im,-0.9458172417006346-0.3246994692046836im,-0.9863613034027224+0.1645945902807336im,-0.9863613034027224-0.1645945902807336im,-1+0im]

# Test that our answers are close to MATLAB's and at least as accurate
zpkfilter_eq(Butterworth(19), ZeroPoleGain(Union{}[], matlab_p, 1))
zpkfilter_accuracy(Butterworth(19), ZeroPoleGain(Union{}[], matlab_p, 1), Butterworth(BigFloat, 19))

#
# Chebyshev type I filter
#

# Poles of 20 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
#=
    [z, p, k] = cheb1ap(20, 1)
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
=#
matlab_p = [-0.005606643513655412+0.9994594480354867im,-0.01668187637027696+0.9748494394091176im,-0.02734634606668882+0.9262354022446964im,-0.03733745796561948+0.854814375477951im,-0.04640919773351807+0.7623449818206758im,-0.05433818902992066+0.6511041246330369im,-0.06092919377429644+0.5238309230062319im,-0.06601991955554824+0.3836592655598524im,-0.06948501580979141+0.2340406437032877im,-0.07123916036725134+0.07865916446238384im,-0.07123916036725134-0.07865916446238384im,-0.06948501580979141-0.2340406437032877im,-0.06601991955554824-0.3836592655598524im,-0.06092919377429644-0.5238309230062319im,-0.05433818902992066-0.6511041246330369im,-0.04640919773351807-0.7623449818206758im,-0.03733745796561948-0.854814375477951im,-0.02734634606668882-0.9262354022446964im,-0.01668187637027696-0.9748494394091176im,-0.005606643513655412-0.9994594480354867im]
matlab_k = 3.748372513504540e-06

# Test that our answers are close to MATLAB's and at least as accurate
zpkfilter_eq(Chebyshev1(20, 1), ZeroPoleGain(Union{}[], matlab_p, matlab_k))
zpkfilter_accuracy(Chebyshev1(20, 1), ZeroPoleGain(Union{}[], matlab_p, matlab_k), Chebyshev1(BigFloat, 20, 1))

# Poles of 19 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
#=
    [z, p, k] = cheb1ap(19, 1)
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
matlab_p = [-0.006212227114275604+0.9994004289494353im,-0.0184672279812169+0.9721393904901743im,-0.03021849100930104+0.9183609236365976im,-0.04114547237637834+0.8395319647744703im,-0.05095011251526942+0.7378027624097928im,-0.0593649664017747+0.6159482239948341im,-0.06616049875775722+0.4772922236864799im,-0.07115134517517649+0.3256169357239265im,-0.07420136837393095+0.1650596665748249im,-0.07522737167197566+0im,-0.07420136837393095-0.1650596665748249im,-0.07115134517517649-0.3256169357239265im,-0.06616049875775722-0.4772922236864799im,-0.0593649664017747-0.6159482239948341im,-0.05095011251526942-0.7378027624097928im,-0.04114547237637834-0.8395319647744703im,-0.03021849100930104-0.9183609236365976im,-0.0184672279812169-0.9721393904901743im,-0.006212227114275604-0.9994004289494353im]
matlab_k = 7.496745027009062e-06

# Test that our answers are close to MATLAB's and at least as accurate
zpkfilter_eq(Chebyshev1(19, 1), ZeroPoleGain(Union{}[], matlab_p, matlab_k))
zpkfilter_accuracy(Chebyshev1(19, 1), ZeroPoleGain(Union{}[], matlab_p, matlab_k), Chebyshev1(BigFloat, 19, 1))

#
# Chebyshev type II filter
#

# 20 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
#=
    [z, p, k] = cheb2ap(20, 1)
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
matlab_z = [0+1.003092198482826im,0-1.003092198482826im,0+1.028415193665209im,0-1.028415193665209im,0+1.082392200292394im,0-1.082392200292394im,0+1.172827696614009im,0-1.172827696614009im,0+1.315086999890785im,0-1.315086999890785im,0+1.539769043222366im,0-1.539769043222366im,0+1.913880855430943im,0-1.913880855430943im,0+2.613125929752753im,0-2.613125929752753im,0+4.283657569731185im,0-4.283657569731185im,0+12.74549484318238im,0-12.74549484318238im]
matlab_p = [-0.001929675700544753-1.002788598382701im,-0.006034873757627809-1.028072311003105im,-0.010957850142167-1.081957626858823im,-0.01756373888305481-1.172213900920771im,-0.02744258171352171-1.314120759169777im,-0.04403133250736632-1.538048178699878im,-0.07621942506983267-1.910267531994771im,-0.1536690174908018-2.6032737531287im,-0.4316587213861887-4.238414904224276im,-3.610083779585366-11.62012078608171im,-3.610083779585366+11.62012078608171im,-0.4316587213861887+4.238414904224276im,-0.1536690174908018+2.6032737531287im,-0.07621942506983267+1.910267531994771im,-0.04403133250736632+1.538048178699878im,-0.02744258171352171+1.314120759169777im,-0.01756373888305481+1.172213900920771im,-0.010957850142167+1.081957626858823im,-0.006034873757627809+1.028072311003105im,-0.001929675700544753+1.002788598382701im]
matlab_k = 0.8912509381337452
matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

# The gain shows a greater discrepancy with the gain of the enhanced
# precision filter than MATLAB's, but AFAICT our filter is still more
# accurate: the poles and zeros are each more accurate, and computing
# the gain in enhanced precision yields the same value.

f = Chebyshev2(20, 1)
zpkfilter_eq(f, matlab_f)
zpkfilter_accuracy(f, matlab_f, Chebyshev2(BigFloat, 20, 1), relerr=4)

# 19 pole Butterworth filter prototype with 1 dB passband ripple from MATLAB 2013b:
#=
    [z, p, k] = cheb2ap(19, 1)
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
matlab_z = [0+1.003427212662145im,0-1.003427212662145im,0+1.031565634068626im,0-1.031565634068626im,0+1.091973276457601im,0-1.091973276457601im,0+1.194505544554793im,0-1.194505544554793im,0+1.359205519207709im,0-1.359205519207709im,0+1.628100459889458im,0-1.628100459889458im,0+2.101072544213108im,0-2.101072544213108im,0+3.079770972368366im,0-3.079770972368366im,0+6.075533820974263im,0-6.075533820974263im]
matlab_p = [-0.002139218568677465-1.003090263923361im,-0.006720708432694556-1.031180124433753im,-0.01232195186345456-1.091472452757043im,-0.02007306444044342-1.193772338103026im,-0.03217434756636293-1.357992987527419im,-0.05375960613144559-1.625783486413104im,-0.09966387580015423-2.09563676642751im,-0.2295218019875045-3.061543704346429im,-0.9149742173123645-5.932401750359961im,-38.84170469908301+0im,-0.9149742173123645+5.932401750359961im,-0.2295218019875045+3.061543704346429im,-0.09966387580015423+2.09563676642751im,-0.05375960613144559+1.625783486413104im,-0.03217434756636293+1.357992987527419im,-0.02007306444044342+1.193772338103026im,-0.01232195186345456+1.091472452757043im,-0.006720708432694556+1.031180124433753im,-0.002139218568677465+1.003090263923361im]
matlab_k = 37.33930783884512
matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

f = Chebyshev2(19, 1)
zpkfilter_eq(f, matlab_f)
zpkfilter_accuracy(f, matlab_f, Chebyshev2(BigFloat, 19, 1), relerr=2)

#
# Elliptic filter
#

# 20 pole elliptic filter prototype with 0.1 dB passband ripple and 10
# dB stopband ripple from MATLAB 2013b:
#=
    [z, p, k] = ellipap(20, 0.1, 10)
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
matlab_z = [0-1.953252853757711im,0+1.953252853757711im,0-1.069599937693626im,0+1.069599937693626im,0-1.007032209276402im,0+1.007032209276402im,0-1.000730367252677im,0+1.000730367252677im,0-1.000076070678288im,0+1.000076070678288im,0-1.000007925871874im,0+1.000007925871874im,0-1.000000826312751im,0+1.000000826312751im,0-1.000000086632384im,0+1.000000086632384im,0-1.000000009576065im,0+1.000000009576065im,0-1.000000001633984im,0+1.000000001633984im]
matlab_p = [-0.6665552331000151-0.9586923119289487im,-0.6665552331000151+0.9586923119289487im,-0.06598422536270995-1.015686365686964im,-0.06598422536270995+1.015686365686964im,-0.006773291804375061-1.001821704927991im,-0.006773291804375061+1.001821704927991im,-0.0007045345890457818-1.000191780070487im,-0.0007045345890457818+1.000191780070487im,-7.339097059946433e-05-1.000020003037497im,-7.339097059946433e-05+1.000020003037497im,-7.646273200288023e-06-1.000002084836051im,-7.646273200288023e-06+1.000002084836051im,-7.966432964216136e-07-1.000000217756124im,-7.966432964216136e-07+1.000000217756124im,-8.299912436196738e-08-1.000000023227775im,-8.299912436196738e-08+1.000000023227775im,-8.637784231770464e-09-1.000000002962719im,-8.637784231770464e-09+1.000000002962719im,-8.070974080610083e-10-1.000000000874029im,-8.070974080610083e-10+1.000000000874029im]
matlab_k = 0.3162277662398871
matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

f = Elliptic(20, 0.1, 10)
zpkfilter_eq(f, matlab_f)
zpkfilter_accuracy(f, matlab_f, Elliptic(BigFloat, 20, 0.1, 10), eps=1e-9)

# 19 pole elliptic filter prototype with 0.1 dB passband ripple and 10
# dB stopband ripple from MATLAB 2013b:
#=
    [z, p, k] = ellipap(19, 0.1, 10)
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
matlab_z = [0-1.232609672486912im,0+1.232609672486912im,0-1.021948255265862im,0+1.021948255265862im,0-1.002264470627064im,0+1.002264470627064im,0-1.000235691798532im,0+1.000235691798532im,0-1.000024555183434im,0+1.000024555183434im,0-1.000002559984621im,0+1.000002559984621im,0-1.000000268394018im,0+1.000000268394018im,0-1.00000002966741im,0+1.00000002966741im,0-1.000000005062212im,0+1.000000005062212im]
matlab_p = [-0.2102804655411482-1.034611529620449im,-0.2102804655411482+1.034611529620449im,-0.0210635609742992-1.005504639644632im,-0.0210635609742992+1.005504639644632im,-0.002183581313750937-1.00059265523267im,-0.002183581313750937+1.00059265523267im,-0.0002273805846473938-1.000061954788962im,-0.0002273805846473938+1.000061954788962im,-2.368886404726722e-05-1.000006458811389im,-2.368886404726722e-05+1.000006458811389im,-2.468065209676509e-06-1.00000067462382im,-2.468065209676509e-06+1.00000067462382im,-2.571378798158274e-07-1.000000071961468im,-2.571378798158274e-07+1.000000071961468im,-2.676054057433227e-08-1.000000009178736im,-2.676054057433227e-08+1.000000009178736im,-2.500451772444974e-09-1.00000000270781im,-2.500451772444974e-09+1.00000000270781im,-1.309071549907812+0im]
matlab_k = 0.9266824319626478
matlab_f = ZeroPoleGain(matlab_z, matlab_p, matlab_k)

f = Elliptic(19, 0.1, 10)
zpkfilter_eq(f, matlab_f)
zpkfilter_accuracy(f, matlab_f, Elliptic(BigFloat, 19, 0.1, 10), eps=4e-9)

#
# Conversion between zpk and tf
#

# Output of [z, p, k] = buttap(20); [b, a] = zp2tf(z, p, k)
m_a = [1,12.74549484318237,81.22381939879423,343.6513712403923,1081.352361133001,2687.409807920676,5468.931438945091,9326.061201886809,13528.36656744904,16852.27707949905,18122.54155403868,16852.27707949905,13528.36656744904,9326.061201886809,5468.931438945092,2687.409807920676,1081.352361133001,343.6513712403923,81.22381939879423,12.74549484318237,1]

f = convert(PolynomialRatio, Butterworth(20))
@test coefb(f) ≈ [1]
@test coefa(f) ≈ m_a

# Test that our answers are more accurate than MATLAB's
accurate_a = coefa(convert(PolynomialRatio, Butterworth(BigFloat, 20)))
@test coefa(f) ≈ map(Float64, accurate_a)
@test sum(abs.(coefa(f) - accurate_a)) <= sum(abs.(m_a - accurate_a))

#
# Conversion between tf and zpk
#

f = convert(ZeroPoleGain, convert(PolynomialRatio, Butterworth(20)))
zpkfilter_eq(f, butter, 1e-6)

# TODO: Make this more accurate!

# Output of [z, p, k] = buttap(20); [b, a] = zp2tf(z, p, k); tf2zpk(b, a)
m_p = [-0.07845909573254482+0.9969173337335029im,-0.07845909573254482-0.9969173337335029im,-0.2334453637958131+0.9723699203918822im,-0.2334453637958131-0.9723699203918822im,-0.3826834327796701+0.9238795325396184im,-0.3826834327796701-0.9238795325396184im,-0.5224985628488221+0.8526401643454914im,-0.5224985628488221-0.8526401643454914im,-0.6494480541398985+0.7604059651905597im,-0.6494480541398985-0.7604059651905597im,-0.760405952587916+0.6494480502272874im,-0.760405952587916-0.6494480502272874im,-0.8526401859847815+0.5224985598169277im,-0.8526401859847815-0.5224985598169277im,-0.9238795057196649+0.3826834415328767im,-0.9238795057196649-0.3826834415328767im,-0.9969173244841298+0.07845911266921719im,-0.9969173244841298-0.07845911266921719im,-0.97236994351794+0.2334453500964366im,-0.97236994351794-0.2334453500964366im]
# println(complex128([sort(f.p, lt=lt) - sort(Butterworth(BigFloat, 20).p, lt=lt) sort(m_p, lt=lt) - sort(Butterworth(BigFloat, 20).p, lt=lt)]))
zpkfilter_accuracy(f, ZeroPoleGain(Union{}[], m_p, 1), Butterworth(BigFloat, 20); eps=1e-6, relerr=9)

#
# Frequency scaling
#

# Output of:
#=
    [z, p, k] = butter(20, 0.5, 's')
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
m_z = Union{}[]
m_p = [-0.498458666866501+0.03922954786337574im,-0.498458666866501-0.03922954786337574im,-0.4861849601989526+0.1167226819283556im,-0.4861849601989526-0.1167226819283556im,-0.4619397662555457+0.191341716182291im,-0.4619397662555457-0.191341716182291im,-0.4263200821771486+0.2612492823580798im,-0.4263200821771486-0.2612492823580798im,-0.3802029827999479+0.3247240241651102im,-0.3802029827999479-0.3247240241651102im,-0.3247240241650893+0.3802029827999696im,-0.3247240241650893-0.3802029827999696im,-0.2612492823580033+0.4263200821770542im,-0.2612492823580033-0.4263200821770542im,-0.1913417161825348+0.461939766255657im,-0.1913417161825348-0.461939766255657im,-0.1167226819279489+0.486184960198832im,-0.1167226819279489-0.486184960198832im,-0.03922954786392431+0.4984586668665641im,-0.03922954786392431-0.4984586668665641im]
m_k = 9.5367431640625e-07
m_f = ZeroPoleGain(m_z, m_p, m_k)

f = analogfilter(Lowpass(0.5), Butterworth(20))
zpkfilter_eq(f, m_f)

# Test that our answers are more accurate than MATLAB's
accurate_ft = analogfilter(Lowpass(0.5), Butterworth(BigFloat, 20))
zpkfilter_eq(f, accurate_ft)
zpkfilter_accuracy(f, m_f, accurate_ft)

# Output of
#=
    [z, p, k] = ellip(20, 0.5, 20, 0.3, 's')
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
m_z = [0+0.6269347368248613im,0-0.6269347368248613im,3.244237104018702e-18+0.3275342448365736im,3.244237104018702e-18-0.3275342448365736im,-6.938893903907228e-18+0.3032929712822612im,-6.938893903907228e-18-0.3032929712822612im,-1.387778780781446e-17+0.3004076906327037im,-1.387778780781446e-17-0.3004076906327037im,-4.163336342344337e-17+0.3000506878920809im,-4.163336342344337e-17-0.3000506878920809im,6.938893903907228e-18+0.3000063060899686im,6.938893903907228e-18-0.3000063060899686im,-6.938893903907228e-18+0.3000007854095975im,-6.938893903907228e-18-0.3000007854095975im,-1.387778780781446e-17+0.30000009864114im,-1.387778780781446e-17-0.30000009864114im,1.387778780781446e-17+0.3000000132267227im,1.387778780781446e-17-0.3000000132267227im,2.775557561562891e-17+0.3000000027630185im,2.775557561562891e-17-0.3000000027630185im]
m_p = [-0.1398677695014979+0.188138982052175im,-0.1398677695014979-0.188138982052175im,-0.02345522399724158+0.2890911154818891im,-0.02345522399724158-0.2890911154818891im,-0.003009381182870904+0.2987256455960404im,-0.003009381182870904-0.2987256455960404im,-0.0003757590363774529+0.2998428279105927im,-0.0003757590363774529-0.2998428279105927im,-4.676623056761031e-05+0.299980469734862im,-4.676623056761031e-05-0.299980469734862im,-5.818096677288886e-06+0.2999975716762768im,-5.818096677288886e-06-0.2999975716762768im,-7.237817530730828e-07+0.2999996988517393im,-7.237817530730828e-07-0.2999996988517393im,-9.003660077222531e-08+0.2999999634684986im,-9.003660077222531e-08-0.2999999634684986im,-1.117941275663448e-08+0.2999999963792452im,-1.117941275663448e-08-0.2999999963792452im,-1.220070558471953e-09+0.300000000410979im,-1.220070558471953e-09-0.300000000410979im]
m_k = 0.1000000000276181
m_f = ZeroPoleGain(m_z, m_p, m_k)

f = convert(ZeroPoleGain, analogfilter(Lowpass(0.3), Elliptic(20, 0.5, 20)))
zpkfilter_eq(f, m_f)

#
# High pass filter construction
#

# Output of
#=
    [z, p, k] = butter(20, 0.5, 'high', 's')
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
m_z = [0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im]
m_p = [-0.4984586668666513+0.03922954786370267im,-0.4984586668666513-0.03922954786370267im,-0.4861849601986916+0.1167226819280185im,-0.4861849601986916-0.1167226819280185im,-0.4619397662557133+0.191341716182565im,-0.4619397662557133-0.191341716182565im,-0.4263200821770351+0.2612492823579622im,-0.4263200821770351-0.2612492823579622im,-0.3802029828000168+0.3247240241650909im,-0.3802029828000168-0.3247240241650909im,-0.3247240241650923+0.3802029828000145im,-0.3247240241650923-0.3802029828000145im,-0.03922954786392298+0.4984586668665641im,-0.03922954786392298-0.4984586668665641im,-0.2612492823579744+0.4263200821770468im,-0.2612492823579744-0.4263200821770468im,-0.1913417161825449+0.4619397662556436im,-0.1913417161825449-0.4619397662556436im,-0.1167226819279526+0.4861849601988384im,-0.1167226819279526-0.4861849601988384im]
m_k = 1
m_f = ZeroPoleGain(m_z, m_p, m_k)

f = analogfilter(Highpass(0.5), Butterworth(20))
zpkfilter_eq(f, m_f)

# Test that our answers are more accurate than MATLAB's
accurate_ft = analogfilter(Highpass(0.5), Butterworth(BigFloat, 20))
zpkfilter_eq(f, accurate_ft)
zpkfilter_accuracy(f, m_f, accurate_ft)

# Output of
#=
    [z, p, k] = ellip(20, 0.5, 20, 0.4, 'high', 's')
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
m_z = [4.289315073341702e-17+0.1914074830304432im,4.289315073341702e-17-0.1914074830304432im,-1.734723475976807e-17+0.3663739040779542im,-1.734723475976807e-17-0.3663739040779542im,-1.387778780781446e-17+0.3956570424057772im,-1.387778780781446e-17-0.3956570424057772im,-2.775557561562891e-17+0.3994571502056489im,-2.775557561562891e-17-0.3994571502056489im,2.775557561562891e-17+0.3999324275609072im,2.775557561562891e-17-0.3999324275609072im,8.673617379884035e-18+0.3999915920567793im,8.673617379884035e-18-0.3999915920567793im,4.163336342344337e-17+0.399998952789945im,4.163336342344337e-17-0.399998952789945im,1.387778780781446e-17+0.3999998684785239im,1.387778780781446e-17-0.3999998684785239im,-2.775557561562891e-17+0.3999999963159753im,-2.775557561562891e-17-0.3999999963159753im,6.938893903907228e-18+0.3999999823643708im,6.938893903907228e-18-0.3999999823643708im]
m_p = [-0.3053922020601836+0.4107892635097733im,-0.3053922020601836-0.4107892635097733im,-0.03345814372376468+0.4123794380387859im,-0.03345814372376468-0.4123794380387859im,-0.004046405082073851+0.4016656239382425im,-0.004046405082073851-0.4016656239382425im,-0.0005015366410947575+0.4002090441157958im,-0.0005015366410947575-0.4002090441157958im,-6.236309209528557e-05+0.400026032326622im,-6.236309209528557e-05-0.400026032326622im,-7.757587819075262e-06+0.4000032376407234im,-7.757587819075262e-06-0.4000032376407234im,-9.650442748948551e-07+0.4000004015290894im,-9.650442748948551e-07-0.4000004015290894im,-1.200488304875513e-07+0.4000000487086385im,-1.200488304875513e-07-0.4000000487086385im,-1.490588337482723e-08+0.4000000048276729im,-1.490588337482723e-08-0.4000000048276729im,-1.626761100825824e-09+0.3999999994520274im,-1.626761100825824e-09-0.3999999994520274im]
m_k = 0.9440608762859236
m_f = ZeroPoleGain(m_z, m_p, m_k)

f = convert(ZeroPoleGain, analogfilter(Highpass(0.4), Elliptic(20, 0.5, 20)))
zpkfilter_eq(f, m_f)

#
# Band pass filter construction
#

# Output of
#=
    [z, p, k] = butter(10, [0.25 0.75], 's')
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
m_z = [0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im,0+0im]
m_p = [-0.05852661385905927+0.7442329445709392im,-0.05852661385905927-0.7442329445709392im,-0.1665549471789192+0.6992510095482473im,-0.1665549471789192-0.6992510095482473im,-0.2479888914449238+0.615606020752747im,-0.2479888914449238-0.615606020752747im,-0.2869902139927978+0.5070583156429217im,-0.2869902139927978-0.5070583156429217im,-0.2738314211701858+0.3979722116506606im,-0.2738314211701858-0.3979722116506606im,-0.2200127491273796+0.3197549791305416im,-0.2200127491273796-0.3197549791305416im,-0.1585130481013887+0.280063065773154im,-0.1585130481013887-0.280063065773154im,-0.1055644991483508+0.2620526301594698im,-0.1055644991483508-0.2620526301594698im,-0.06044030269085397+0.2537477474540631im,-0.06044030269085397-0.2537477474540631im,-0.01969061866105566+0.2503887742733703im,-0.01969061866105566-0.2503887742733703im]
m_k = 0.0009765624999999997
m_f = ZeroPoleGain(m_z, m_p, m_k)

f = analogfilter(Bandpass(0.25, 0.75), Butterworth(10))
zpkfilter_eq(f, m_f)

# Test that our answers are more accurate than MATLAB's
accurate_ft = analogfilter(Bandpass(0.25, 0.75), Butterworth(BigFloat, 10))
zpkfilter_eq(f, accurate_ft)
zpkfilter_accuracy(f, m_f, accurate_ft)

# Output of
#=
    [z, p, k] = ellip(10, 0.5, 20, [0.2 0.8], 's')
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
m_z = [1.387778780781446e-16+1.370743857088904im,1.387778780781446e-16-1.370743857088904im,-5.204170427930421e-18+0.1167249440313361im,-5.204170427930421e-18-0.1167249440313361im,0+0.8445880526167882im,0-0.8445880526167882im,-1.734723475976807e-18+0.8053338099984824im,-1.734723475976807e-18-0.8053338099984824im,-1.387778780781446e-17+0.1894414673571001im,-1.387778780781446e-17-0.1894414673571001im,0+0.1986753790956593im,0-0.1986753790956593im,5.551115123125783e-17+0.800148411646084im,5.551115123125783e-17-0.800148411646084im,1.108054620280186e-16+0.8007109707175198im,1.108054620280186e-16-0.8007109707175198im,8.673617379884035e-18+0.199822415142662im,8.673617379884035e-18-0.199822415142662im,4.85722573273506e-17+0.1999629039703328im,4.85722573273506e-17-0.1999629039703328im]
m_p = [-0.2019763467382846+0.6120628823384283im,-0.2019763467382846-0.6120628823384283im,-0.03721004055421873+0.7823115667177283im,-0.03721004055421873-0.7823115667177283im,-0.004810677763557675+0.7980130970318471im,-0.004810677763557675-0.7980130970318471im,-0.0006000695021327641+0.7998051371015736im,-0.0006000695021327641-0.7998051371015736im,-6.552554548894252e-05+0.8000220660203423im,-6.552554548894252e-05-0.8000220660203423im,-0.07779252812598642+0.235740074311183im,-0.07779252812598642-0.235740074311183im,-0.009705975806501373+0.2040604370920958im,-0.009705975806501373-0.2040604370920958im,-0.001208621817121024+0.2004906765377495im,-0.001208621817121024-0.2004906765377495im,-0.0001500903998086778+0.2000486149853388im,-0.0001500903998086778-0.2000486149853388im,-1.638048261950673e-05+0.1999944823054316im,-1.638048261950673e-05-0.1999944823054316im]
m_k = 0.1000000000000025
m_f = ZeroPoleGain(m_z, m_p, m_k)

f = convert(ZeroPoleGain, analogfilter(Bandpass(0.2, 0.8), Elliptic(10, 0.5, 20)))
zpkfilter_eq(f, m_f)

#
# Band stop filter construction
#

# Output of
#=
    [z, p, k] = butter(10, [0.25 0.75], 'stop', 's')
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
m_z = [0+0.4330127018922193im,-0-0.4330127018922193im,0+0.4330127018922193im,-0-0.4330127018922193im,0+0.4330127018922193im,-0-0.4330127018922193im,0+0.4330127018922193im,-0-0.4330127018922193im,0+0.4330127018922193im,-0-0.4330127018922193im,0+0.4330127018922193im,-0-0.4330127018922193im,0+0.4330127018922193im,-0-0.4330127018922193im,0+0.4330127018922193im,-0-0.4330127018922193im,0+0.4330127018922193im,-0-0.4330127018922193im,0+0.4330127018922193im,-0-0.4330127018922193im]
m_p = [-0.05852661385905998+0.7442329445709392im,-0.05852661385905998-0.7442329445709392im,-0.1665549471789191+0.6992510095482469im,-0.1665549471789191-0.6992510095482469im,-0.2479888914449218+0.6156060207527428im,-0.2479888914449218-0.6156060207527428im,-0.2869902139927988+0.5070583156429307im,-0.2869902139927988-0.5070583156429307im,-0.2738314211701937+0.3979722116506589im,-0.2738314211701937-0.3979722116506589im,-0.2200127491273721+0.3197549791305355im,-0.2200127491273721-0.3197549791305355im,-0.1585130481013898+0.2800630657731571im,-0.1585130481013898-0.2800630657731571im,-0.1055644991483515+0.2620526301594707im,-0.1055644991483515-0.2620526301594707im,-0.06044030269085383+0.2537477474540633im,-0.06044030269085383-0.2537477474540633im,-0.01969061866105538+0.2503887742733704im,-0.01969061866105538-0.2503887742733704im]
m_k = 1
m_f = ZeroPoleGain(m_z, m_p, m_k)

f = analogfilter(Bandstop(0.25, 0.75), Butterworth(10))
zpkfilter_eq(f, m_f)

# Test that our answers are more accurate than MATLAB's
accurate_ft = analogfilter(Bandstop(0.25, 0.75), Butterworth(BigFloat, 10))
zpkfilter_eq(f, accurate_ft)
zpkfilter_accuracy(f, m_f, accurate_ft)

#=
    [z, p, k] = ellip(10, 0.5, 20, [0.2 0.8], 'stop', 's')
    sprintf('%.16g%+.16gim,', [real(z) imag(z)]')
    sprintf('%.16g%+.16gim,', [real(p) imag(p)]')
    sprintf('%.16g', k)
=#
m_z = [-1.402517154063671e-17+0.5685129776423937im,-1.402517154063671e-17-0.5685129776423937im,-6.156235460644255e-18+0.2814359676774932im,-6.156235460644255e-18-0.2814359676774932im,-2.775557561562891e-17+0.760016963542146im,-2.775557561562891e-17-0.760016963542146im,-8.326672684688674e-17+0.7947386855861924im,-8.326672684688674e-17-0.7947386855861924im,5.551115123125783e-17+0.799851645232653im,5.551115123125783e-17-0.799851645232653im,8.326672684688674e-17+0.7992903327101535im,8.326672684688674e-17-0.7992903327101535im,0+0.2105216168522108im,0-0.2105216168522108im,-4.423544863740858e-17+0.2013240363176547im,-4.423544863740858e-17-0.2013240363176547im,6.938893903907228e-18+0.200037095570968im,6.938893903907228e-18-0.200037095570968im,5.551115123125783e-17+0.2001775743458425im,5.551115123125783e-17-0.2001775743458425im]
m_p = [-0.3778783653696471+0.7819846151027801im,-0.3778783653696471-0.7819846151027801im,-0.04044889726976308+0.8144610469677399im,-0.04044889726976308-0.8144610469677399im,-0.004859687206700769+0.8019364482578354im,-0.004859687206700769-0.8019364482578354im,-0.0006006726587018951+0.8001940294215477im,-0.0006006726587018951-0.8001940294215477im,-6.551807382872143e-05+0.799977924147539im,-6.551807382872143e-05-0.799977924147539im,-0.08015534922525629+0.1658741427311564im,-0.08015534922525629-0.1658741427311564im,-0.009732315617052817+0.1959655887284332im,-0.009732315617052817-0.1959655887284332im,-0.00120901710644529+0.1995097303565812im,-0.00120901710644529-0.1995097303565812im,-0.0001500952639848091+0.1999513917362522im,-0.0001500952639848091-0.1999513917362522im,-1.638042236460496e-05+0.2000055177738602im,-1.638042236460496e-05-0.2000055177738602im]
m_k = 0.9440608762859234
m_f = ZeroPoleGain(m_z, m_p, m_k)

f = analogfilter(Bandstop(0.2, 0.8), Elliptic(10, 0.5, 20))
zpkfilter_eq(f, m_f)

#
# IIR digital filter creation
#

# Output of [b, a] = butter(20, 0.5)
m_b = [8.555155699386467e-06,0.0001711031139877293,0.001625479582883429,0.009752877497300572,0.04144972936352743,0.1326391339632878,0.3315978349082195,0.663195669816439,1.077692963451713,1.436923951268951,1.580616346395846,1.436923951268951,1.077692963451713,0.663195669816439,0.3315978349082195,0.1326391339632878,0.04144972936352743,0.009752877497300572,0.001625479582883429,0.0001711031139877293,8.555155699386467e-06]
m_a = [1,-9.216626485446244e-15,2.719352210996272,-2.276898740746355e-14,2.954502383565172,-2.181204855972978e-14,1.663404391798304,-9.886771142150971e-15,0.5272563634534938,-1.578399203426284e-15,0.09588816575455775,3.201156931050434e-16,0.009788250615318749,1.651392612127069e-16,0.0005259945756568433,3.034673947897289e-17,1.306795536212839e-05,3.205637858266449e-18,1.137794123801438e-07,9.233765069631315e-20,1.463813780864114e-10]
m_f = PolynomialRatio(m_b, m_a)

f = convert(PolynomialRatio, digitalfilter(Lowpass(0.5), Butterworth(20)))
tffilter_eq(f, m_f)

# Output of [b, a] = butter(20, 0.3, 'high')
m_b = [0.001797932208654781,-0.03595864417309561,0.3416071196444083,-2.04964271786645,8.710981550932413,-27.87514096298372,69.6878524074593,-139.3757048149186,226.4855203242427,-301.980693765657,332.1787631422226,-301.980693765657,226.4855203242427,-139.3757048149186,69.6878524074593,-27.87514096298372,8.710981550932413,-2.04964271786645,0.3416071196444083,-0.03595864417309561,0.001797932208654781]
m_a = [1,-7.99378842342338,32.24776489137069,-86.15266953725099,169.613151737271,-260.1393113842562,321.1838937487486,-325.8437287835835,275.1584488419961,-194.8984651152305,116.2180984348268,-58.35420576527059,24.59604119108612,-8.645276014889031,2.506713994242244,-0.5898102957633149,0.1098868139124379,-0.01561185396257133,0.001590100744530764,-0.0001034620065564712,3.232560226924501e-06]
m_f = PolynomialRatio(m_b, m_a)

f = convert(PolynomialRatio, digitalfilter(Highpass(0.3), Butterworth(20)))
tffilter_eq(f, m_f)

# Output of [b, a] = butter(7, [0.1 0.4])
m_b = [0.0009628940476886417,0,-0.006740258333820491,0,0.02022077500146148,0,-0.03370129166910246,0,0.03370129166910246,0,-0.02022077500146148,0,0.006740258333820491,0,-0.0009628940476886417]
m_a = [1,-7.763446908073567,29.02150968091795,-69.51575164764138,119.5092528789707,-156.1341322158939,159.8101214681526,-130.0847734566226,84.56503905088692,-43.67462798637182,17.64013826588496,-5.405682857208016,1.18954580728982,-0.1686342472129378,0.01166272804918956]
m_f = PolynomialRatio(m_b, m_a)

f = convert(PolynomialRatio, digitalfilter(Bandpass(0.1, 0.4), Butterworth(7)))
tffilter_eq(f, m_f)

# Output of [b, a] = butter(4, [0.2 0.25], 'stop')
m_b = [0.8142545568862771,-4.968628804996451,14.62659178092118,-26.46885325790832,32.03454308505385,-26.46885325790833,14.62659178092118,-4.968628804996451,0.8142545568862771]
m_a = [1,-5.789124865012897,16.16849931613174,-27.77244867328942,31.92014561624699,-25.06018377653858,13.16458034390183,-4.253206810966346,0.6630104843858931]
m_f = PolynomialRatio(m_b, m_a)

f = convert(PolynomialRatio, digitalfilter(Bandstop(0.2, 0.25), Butterworth(4)))
tffilter_eq(f, m_f)

# Output of [b, a] = ellip(10, 0.7, 13, 0.4)
m_b = [0.2797942479261318,-0.546743606764335,1.704005589668604,-2.159911306584612,3.642222926761125,-3.210370872716824,3.642222926761128,-2.159911306584614,1.704005589668607,-0.5467436067643362,0.2797942479261327]
m_a = [1,-3.029857317356546,7.982710263075236,-12.64912738840943,17.41532979440429,-17.06917827530357,14.5191214973975,-8.798609366466101,4.541539167299845,-1.419553904214128,0.3565809079146795]
m_f = PolynomialRatio(m_b, m_a)

f = convert(PolynomialRatio, digitalfilter(Lowpass(0.4), Elliptic(10, 0.7, 13)))
tffilter_eq(f, m_f)

# Output of [b, a] = ellip(11, 0.7, 17, 0.2, 'high')
m_b = [0.5737867105773529,-5.288884465999311,23.07976318423696,-62.6831981311594,117.4527496444787,-159.209133733912,159.2091337339121,-117.4527496444788,62.68319813115945,-23.079763184237,5.288884465999319,-0.5737867105773541]
m_a = [1,-8.215565996888206,32.26875224082856,-79.26375768096621,134.6328215735552,-165.4634116073934,149.7532535390064,-99.57881504158784,47.56023313164873,-15.48547977535943,3.075141867250689,-0.2777992862430019]
m_f = PolynomialRatio(m_b, m_a)

f = convert(PolynomialRatio, digitalfilter(Highpass(0.2), Elliptic(11, 0.7, 17)))
tffilter_eq(f, m_f)

# Output of [b, a] = ellip(10, 0.5, 2, [0.2 0.3])
m_b = [0.7571867287399742,-10.55194327196835,73.54073467886568,-338.6101060102045,1150.502868412841,-3057.625219554799,6581.673203566225,-11732.90330745444,17574.06211811366,-22320.04109078996,24158.39342879073,-22320.04109078995,17574.06211811365,-11732.90330745443,6581.673203566216,-3057.625219554793,1150.502868412839,-338.6101060102037,73.54073467886548,-10.55194327196831,0.7571867287399715]
m_a = [1,-13.88122062422027,96.34606972585827,-441.7050045576331,1494.048255114342,-3952.085137071483,8465.769048060451,-15015.81246941687,22374.64201809294,-28265.00095948501,30424.52053126076,-27950.27646663918,21879.36223420588,-14520.36940748522,8095.785094725445,-3737.670491911818,1397.479471832041,-408.6480215434079,88.17084042962101,-12.56723162073412,0.8957646298422474]
m_f = PolynomialRatio(m_b, m_a)

f = convert(PolynomialRatio, digitalfilter(Bandpass(0.2, 0.3), Elliptic(10, 0.5, 2)))
tffilter_eq(f, m_f)

# Output of [b, a] = ellip(15, 1.8, 52, [0.6 0.7], 'stop')
m_b = [0.4754470525338341,6.428928219976483,47.47972043083799,245.8985324356427,989.107195159456,3259.995646867947,9099.649115772192,21987.91200123768,46709.48715187636,88217.96445932264,149368.0297431405,228141.9056208539,315792.8309895032,397458.2884412499,455878.980369108,477133.9777299722,455878.9803691082,397458.2884412502,315792.8309895035,228141.9056208542,149368.0297431408,88217.96445932281,46709.48715187647,21987.91200123774,9099.649115772219,3259.995646867959,989.1071951594597,245.8985324356436,47.4797204308382,6.428928219976513,0.475447052533836]
m_a = [1,12.89634891736776,90.77856214451822,447.9765044361528,1716.574771803699,5388.380608104304,14321.19101752449,32940.21064000877,66586.59389679749,119619.3933180238,192555.4138276511,279453.5193544072,367301.6888247615,438620.3975638363,476900.358709065,472649.3359368506,427104.9678341324,351676.250633706,263451.4609235193,179107.598840823,110104.024869051,60897.74794125407,30103.22741782092,13182.01328029648,5052.946545447418,1668.118080642225,463.5003850459777,104.731419523004,18.20949251383043,2.194519514615899,0.1423868039113367]
m_f = PolynomialRatio(m_b, m_a)

f = convert(PolynomialRatio, digitalfilter(Bandstop(0.6, 0.7), Elliptic(15, 1.8, 52)))
tffilter_eq(f, m_f)

#
# Special filters
#

# Output of [num, den] = iirnotch(0.1, 0.1)
notch = iirnotch(0.1, 0.1)
@test notch.b0 ≈ 0.863271264002681
@test notch.b1 ≈ -1.642039521920206
@test notch.b2 ≈ 0.863271264002681
@test notch.a1 ≈ -1.642039521920206
@test notch.a2 ≈ 0.726542528005361

notch = iirnotch(60, 1, fs=1000)
@test notch.b0 ≈ 0.996868235770807
@test notch.b1 ≈ -1.853729290297204
@test notch.b2 ≈ 0.996868235770807
@test notch.a1 ≈ -1.853729290297204
@test notch.a2 ≈ 0.993736471541615

#
# IIR digital filter error conditions
#

# zero pole
@test_throws ErrorException Butterworth(0)
@test_throws ErrorException Chebyshev1(0, 1)
@test_throws ErrorException Chebyshev2(0, 1)
@test_throws ErrorException Elliptic(0, 0.1, 0.5)

# too much passband ripple
@test_throws ErrorException Elliptic(2, 0.5, 0.1)

# negative ripple
@test_throws ErrorException Chebyshev1(2, -1)
@test_throws ErrorException Chebyshev2(2, -1)
@test_throws ErrorException Elliptic(2, 0.1, -0.5)

#
# Specification of sample rate
#

@test Lowpass(20, fs=100) == Lowpass(0.4)
@test Highpass(20, fs=100) == Highpass(0.4)
@test Bandpass(20, 21, fs=100) == Bandpass(0.4, 0.42)
@test Bandstop(20, 21, fs=100) == Bandstop(0.4, 0.42)

#
# kaiserord
#

kaiserord_jl    = kaiserord(0.1, 60)
kaiserord_scipy = (74, 1.7994885471673767)    # (n,beta)=kaiserord(60,.1); (n,beta/pi)
@test kaiserord_jl[1] == kaiserord_scipy[1]
@test kaiserord_jl[2] ≈ kaiserord_scipy[2]

kaiserord_jl    = kaiserord(0.1, 21.0 + eps(21.0))
kaiserord_scipy = (20, 3.087704731700782e-07) # (n,beta)=kaiserord(21.000000000000004,.1); (n,beta/pi)
@test kaiserord_jl[1] == kaiserord_scipy[1]
@test kaiserord_jl[2] ≈ kaiserord_scipy[2]

kaiserord_jl    = kaiserord(0.1, 15)
kaiserord_scipy = (11, 0.0)                   # (n,beta)=kaiserord(15,.1); (n,beta/pi)
@test kaiserord_jl[1] == kaiserord_scipy[1]
@test kaiserord_jl[2] ≈ kaiserord_scipy[2]

#
# Window FIR filter taps
#

winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(128), scale=false))
# firwin(128, 0.25, nyq=.5, scale=False)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_lowpass_fc0.25_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ winfirtaps_scipy

winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(129), scale=false))
# firwin(129, 0.25, nyq=.5, scale=False)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_lowpass_fc0.25_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ winfirtaps_scipy

@test_throws ArgumentError digitalfilter(Highpass(0.25; fs=1),FIRWindow(hamming(128), scale=false))

winfirtaps_jl    = digitalfilter(Highpass(0.25; fs=1),FIRWindow(hamming(129), scale=false))
# firwin(129, 0.25, nyq=.5, scale=False, pass_zero=False)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_highpass_fc0.25_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ winfirtaps_scipy

winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(128), scale=false))
# firwin(128, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_bandpass_fc0.1_0.2_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ winfirtaps_scipy

winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=false))
# firwin(129, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandpass_fc0.1_0.2_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ winfirtaps_scipy

@test_throws ArgumentError digitalfilter(Bandstop(0.1, 0.2; fs=1),FIRWindow(hamming(128), scale=false))

winfirtaps_jl    = digitalfilter(Bandstop(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=false))
# firwin(129, [0.1, 0.2], nyq=.5, scale=False)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandstop_fc0.1_0.2_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ vec(winfirtaps_scipy)

winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(128), scale=true))
# firwin(128, 0.25, nyq=.5)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_lowpass_scaled_fc0.25_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ winfirtaps_scipy

winfirtaps_jl    = digitalfilter(Lowpass(0.25; fs=1),FIRWindow(hamming(129), scale=true))
# firwin(129, 0.25, nyq=.5)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_lowpass_scaled_fc0.25_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ winfirtaps_scipy

winfirtaps_jl    = digitalfilter(Highpass(0.25; fs=1),FIRWindow(hamming(129), scale=true))
# firwin(129, 0.25, nyq=.5, pass_zero=False)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_highpass_scaled_fc0.25_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ winfirtaps_scipy

winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(128), scale=true))
# firwin(128, [0.1, 0.2], nyq=.5, pass_zero=False)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_bandpass_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ winfirtaps_scipy

winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=true))
# firwin(129, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandpass_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ winfirtaps_scipy

winfirtaps_jl    = digitalfilter(Bandstop(0.1, 0.2; fs=1),FIRWindow(hamming(129), scale=true))
# firwin(129, [0.1, 0.2], nyq=.5, scale=False)
winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandstop_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
@test winfirtaps_jl ≈ vec(winfirtaps_scipy)
